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Abstract 

Some results of test runs on a 6^ x 12 lattice with Wilson quarks and gauge group 
SU(2) for a previously proposed fermion algorithm by A. Slavnov are presented. 

1 Introduction. 

The incorporation of fermionic degrees of freedom in lattice simulations presents a serious 
difficulties which stem from their being anticommuting variables. For the most actions 
in use which are quadratic in fermionic fields, one can eliminate fermions by an analytic 
integration. Unfortunately, the resulting expressions involve the determinant of a very 
large matrix, which induces nonlocal long-range interactions on the gauge fields, making 
the numerical simulations rather expensive. 

The most often used method to simulate fermionic field theories is now Hybrid Monte 
Carlo (HMC) algorithm which addresses the problem of nonlocality by linearizing the 
action in a succession of molecular dynamics steps. It is exact and easily implementable, 
but still the fermionic simulations turn out to be very time consuming and alternative 
approaches are welcome. 

An interesting idea has been put forward by M. Liischer 0, who proposed to ap- 
proximate the inverse of the fermion matrix with a polynomial, and then interprete the 
determinant as the partition function of a local bosonic system. This technique has been 
studied thorougly for the last five years and it was claimed to be competitive to the HMC 
algorithm and even superior when fermions are heavy (a detailed analysis can be found in 
Refspl, However the efficiency of Liiescher's method is abated when quark mass 



decreases. 

Another approach to the problem of dynamical fermions simulation was proposed by 
A. Slavnov in papers |^, where a D-dimensional fermion determinant was presented 
as a path integral of a (D-l-l)-dimensional local bosonic action. In Ref.0 this procedure 
was tested by numerical simulation of a one dimensional toy model and it was shown that 
correct and accurate results can be obtained with a reasonable size of lattice in auxiliary 
dimension. In this paper we present some results from test runs for SU(2) QCD with 2 
degenerate heavy Wilson quarks, using the version of Slavnov's algorithm from Ref. [0]. 



2 The local bosonic algorithm. 

For the local gauge theory with two flavours of Wilson quarks the effective distribution 
of the gauge field U is given by 

P[U] oc det{D[U] + m)2e-^«[^] (1) 
where D[U] is Wilson-Dirac operator 

^ = ^Ei7.(v, + v;)-v;v,] (2) 

Sg[U] denotes the usual plaquette action and m is the bare quark mass which is related 
to hopping parameter k through k = 1/(8 + 2m). The lattice spacing has been set equal 
to unity for simplicity. It will be advantageous for us to work with hermitean operator 
5 = 75(D + m), so we rewrite expression (|1|) in the form 

P[U] oc detB[U]^e-^^^^^ (3) 

Now following Ref.[^, we introduce five dimensional bosonic fields 4>{x,t) = 0„(x), 
where the extra coordinate t is defined on one dimensional chain of the length L with the 
lattice spacing b 

t = nb ; 0<n<N ; L = Nb (4) 

and four dimensional bosonic fields x{^)- After that we present the determinant in eq.(^ 
in the following form 

detB\Uf = lim det{B[Uf + n^) = lim lim / e--'^'^'>''^^^''^'^^D(f)* D(j)DxDx (5) 
where S^^b,L is a local bosonic action 

N-l 
X n=l 

+b{t<j)l^,{x)BMx) + h.c) + b^<f)l{x)B^Mx) + 

+bexp{~iJ,bn}{x*ix){iJ, + tB)(f)n{x) + h.c.) + 7^ E (6) 
and the free boundary conditions in t for the fields (pn are imposed 

00 = ; 0jv = O (7) 

Let us prove eq.(^), which shows that lattice QCD is a limit of purely bosonic theory 
(see Ref. 0). The bosonic action (^) is a linearized version of the expression in the 
exponent of the following integral 



I,,bAU] = I exp{5: Y: [(C;iexp{-25.6}C + h.c. - 2^*0 - 

a n=l 

-6exp{-/i6n}(x"*(/i + zi?a)C + /^.c)] - (8) 
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where instead of ^-representation we used a basis formed by the eigenvectors of the oper- 
ator B, Ba being corresponding eigenvalues. Indeed, expanding exp{—zBab} in eq.(|^) in 
a Taylor series, keeping only the terms of the order 0{b'^) and substituting 



^''''nllBlrn + h.C.)^b'r^*Blrn (9) 
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we get the expression (|^). The substitution (^) also introduces corrections of the order 
(9(6^) and ensures that the action (H) is strictly positive. 
By changing variables 

0°^exp{-2i?,n6}< ; exp{zi?,n6}< (10) 

we can rewrite the equation as the gaussian integral over with the quadratic form 
which does not depend on B^ 



N-l 



I,,tAU] = f exp{^ [i<PZi€ + h.c. - 20rC) - 

•' a n=l 

-6(exp{-(^ + zB^)bn}x"*{fi + iB^)^^ + h.c)] - ^ E x"*x"}^</>*^</>^X*^X 



Let us integrate the expression (11) over the fields 0. To do so it is sufficient to find 
a stationary point of the exponent. For small b the sum over n can be replaced by the 
integral over continuous variable t = nb and the equations for the stationary point acquire 
the form 



no) = r{L) = (12) 

The replacement of the sum by the integral also introduces corrections of order 0{b^). 
The solutions of these equations are: 

0°(t) = ^^!_^(e-(^-'^<^)* + 1(1 - e-('^-^^-)^) - l) (13) 
b[/i — iBa) ^ L ' 



Substituting these solutions to the integrand ([Tl|), integrating over t and rescaling the 
fields X ~^ VbLx , we get 



limJ^,,,4f/] = / exp{-J2^QQ{l-2e-^'^cosB^L + e-'^'^)}Dx*Dx (14) 
Therefore: 



lim I^,tAU] = det{B[U]^ + /i^) (15) 

6— >0,L— >oo 

The equality (^ is proven. For finite b, L and fi this equation is approximate 
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DetB[U]'^ ^ J e-^^''''^^^''l'^^^D(t)*D(t)Dx*Dx (16) 
and has to be corrected by the terms 

0{b') ; 0(e-^^) ; 0(/^^) (17) 

It was shown that lattice QCD is a hmit of a purely bosonic local theory with the 
action 

S,ff[U, 0, x] = Sg[U] + S^,hAU, 0, X] (18) 

Making use of this action, we can simulate the theory by locally updating the boson fields 
U, (j) and X- 

Choosing the number of points in auxiliary dimension we can fully control the 
systematic errors of the algorithm (p!7|). Indeed, fixing /iL = e we have 

iV=f (19) 

Increasing N we can increase e and decrease h and /x, making the errors (0) small enough 
in comparison with an available level of statistical precision. Since the computational cost 
of the algorithm depends on N , the parameters e, /x, h should be optimally tuned to make 
N as small as possible. 

3 Test results. Investigation of systematic errors of 
the algorithm. 

In this section we present the results from numerical simulations of the bosonic theory ([T8|) 
on a 6^ X 12 lattice with gauge group SU(2) and bare parameters /? = 2.12 and k = 0.15. 
The simulations for the same model and parameters were already done in Ref.p| within 
the framework of Luescher's formulation of the dynamical quarks problem . That will 
allow us to compare the algorithm explored in this paper with the Luescher's one. 

The implementation was done on a Quadrics Q4 machine (APE computer) with 32 
nodes. A full update cycle involved 1 heatbath and 1 overrelaxation sweep for the fields 0„ 
followed by 1 heatbath sweep for the fields x gauge fields U. (Our preparatory tests 
showed that overrelaxing the gauge fields and the fields x does not decrease autocorrelaton 
times Tint- On the contrary, overrelaxing the fields 0„ can decrease r^nt substantially.) For 
each set of parameters /i, b, N we have performed at least 30000 and sometimes up to 
80000 such cycles after thermalization of the system, using the random number generator 
of Ref.||. 

Our test runs were structurally divided into a sequence of "experiments" in which 
the systematic errors (|17D were studied separately. In each "experiment" we fixed some 
of the parameters /iL, holding the corresponding errors fixed and altered the other 
parameters. This way of investigation of the systematic errors allows to determine the 
values of parameters at which one may be able to get correct results at given level of 
statistics, using the bosonic theory with relatively small number of points in auxiliary 
dimension A^. 
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In the first "experiment" we fixed /xL = 4.0, /x = 1.0 to study the alteration of the 
plaquette value when the lattice spacing in auxiliary dimension b changes. By doing so we 
investigated the algorithm's error of the order 0{lP') for fixed errors 0{e~^^) and 
The results of the simulations are reported in Table 1. In Fig. 1 the plaquette values 
multiplied by 100 are plotted against the inverse lattice spacing . The plot shows that 
the plaquette value stabilizes quickly when h decreases. From the data obtained we see 
that at this level of statistics one may be able to do with h = 0.08 to make the systematic 
error of the order 0{b'^) small. 
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20 
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5 
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10 


12.5 


15 




plaquette 


0.5467(2) 


0.5697(3) 


0.5782(4) 


0.5807(8) 


0.5812(12) 


0.5796(6) 




24(2) 


67(8) 


94(18) 


* 


* 





Table 1: Simulation results for the plaquette at /x = 1.0 and fibN = 4.0. The autocorre- 
lation time for the plaquette is measured in unit of iteration sweeps. A star means that the 
autocorrelation time was too long to be measured accurately. In the last column we adduce for 
comparison a value obtained from the Hybrid Monte Carlo for the same lattice and parameters 
/?, k in Ref.§. 



In the second "experiment" we attempted to investigate the character of systematic 
deviations of the order 0{e~^^). For that purpose we fixed = 1 and h = 0.1 in order 
to fix the deviations connected with the errors of the order 0{^'^) and 0{h'^). As a work 
hypothesis we assumed that the correction terms of the order 0{h'^) to the approximation 
(p!6|) do not depend on N at fixed h (it will be seen from the simulation results that this 
hypothesis is wrong). 

The results are reported in Table 2. In Fig. 2 the plaquette values multiplied by 100 
are plotted against the number of points in auxiliary dimension A^. 
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Table 2: Simulation results at/x = 1.0,6 = 0.1. 



The plot shows that the plaquette value does not stabilize when increases. Since the 
error of the order 0(/U^) can not depend on N (the approximation DetB^ ^ Det{B^ + /x^) 
was used even before the bosonization procedure and introduction of auxiliary dimension), 
the only explanation to this result is the dependence on of the error 0{b'^) for fixed b. 

To demonstrate this on a theoretical level, let us return to the expression (^) and 
integrate it over Dcpn for finite b (i.e. not using the simplifying assumption 6 — > 0). Then 
the equations for the stationary point acquire the form 



0n+l 



C-1 + &X°(/i - 25„)e-('^-*^'^)"'' = 

00 = 0Af = 

The solutions of these equations are 



(20) 



211 



6(/i — tBa) ^ ' ' I 2 cosh[/i6 — %bBc\ — 1 

Substituting these solutions to the integrand (|Tl|) , omiting the terms of the order 0{e~^^ 



summing over n, rescaling the fields x ~^ vbLx and keeping only the terms of the order 
0(6^), we get 



I[B]^ /exp{-^ 



1 + ^'^^l + ^(3/i^ - Bl){Bl + f^')] }Dx*Dx (22) 
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From the expression (p2D we see that the algorithm's error terms of the order 0{h'^), 
which are associated to replacement of the sum over n by the integral, explicitly depend 
on N for fixed h. Hence a more careful approach to an analysis of systematic error of the 
order is needed. 

As a rule for interesting cases (and for the model under consideration) the following 
condition holds 

(23) 



Therefore, from the expresions ([T^j2^) one sees that the relative systematic error for 



measurement of arbitrary quantity in Slavnov's algorithm can be represented as follows 



A = Ai + A2 + A3 + A4 

/V7)3 

Ai = ; A2 = G ; A3 = He-^'^ ; A4 = P/x' (24) 

where the values F, G, H, P are approximately constant and depend on the parameters 
b, N, /i only in the next order of perturbation theory. 

Using the formula (0), let us implement a more correct investigation of the algorithm's 
deviation of the order 0{e~^^). 

In the third "experiment" we fixed /i = 1, = 0.04 (A2 = const, A4 = const), and 
increased N (decreasing the errors Ai and A3). The results are reported in Table 3. 
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0.5770(10) 



Table 3: Simulation results at /x = 1.0 , ^ = 0.04. 



From Table 3 one sees that the plaquette value changes slowly when increases 
and it's alteration is not monotonic. Such a behavior can be explained by a mutual 
compensation of the errors Ai and A3. Puting together the results from tables 1 and 3 
one can conclude that at given level of statistics it is enough to hold fiL = 5 in order to 
make the systematic error of the order 0{e~^^) small. 

In the fourth "experiment" we fixed = 0.04, /iL = 3.3 (A2 = const, A3 = const), 
and decreased fi (decreasing the errors Ai and A4). The plaquette value and the masses 
of TT and p mesons were mesured. The results of the simulations are reported in Table 4. 
In Fig. 3 the plaquette values multiplied by 100 are plotted against the value of auxiliary 
"mass" fi. 
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Table 4: Simulation results for the plaquette and meson masses at — — = 0.04, fibN = 3.3 . 



The simulation results for the plaquette and meson masses from Table 4 show that 
the measured values stabilize when /i decreases. From the data obtained we see that at 
this level of statistics one may be able to do with fi = 0.8 to make the error of the order 
0(/i2) small. 

Let us make a conclusion to this section. We investigated the systematic errors of 
Slavnov's algorithm and got the values of algorithm's parameters at which in the tested 
model the errors can be considered small at the given level of statistical precision. These 
values are as follows: 

< 0.03 ; b < 0.08 ; fibN > 5.0 ; fi < 0.8 (25) 



From the conditions (^) we infer that it is enough to take N = 100 points in auxiliary 
dimension to get correct results in the tested model. To check up this statement we made 
a control test run at = 100 ; /i = 0.8 ; b = 0.062. The measured plaquette value 
< P >= 0.5805(11) within the statistical error coincided with the result obtained from 
HMC in paper [|: < P >hmc= 0.5796(6). 
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4 Conclusion. 



In this paper we performed the first investigation of systematic error effects of Slavnov's 
algorithm in a reahstic model. It was shown that at this level of statistics one may be able 
to do with = 100 points in auxiliary dimension. This is comparable with the number 
of bosonic fields which was necessary to get the correct results in the same model from 
Liischer's algorithm in Ref. 0. 

Note that the action (|^) and the effective bosonic action in Liischer's formulation 
of the dynamical fermions problem are similar in a sense that computational effort per 
update cycle is almost the same for both algorithms if the number of bosonic fields in 
Liischer's algorithm is equal to the number of points in auxiliary dimension in action (^) 
(the update of the fields x iii the action @ takes less than 1% of total computer time). 
Accordingly both algorithms are comparable in an amount of computer time which is 
necessary per full update cycle. 

We observed the growth of autocorrelation times when N increases. It seems that 
such an autocorrelation behavior is the general feature of multibosonic algorithms (the 
theoretical justification of this fact can be found in Ref.[0). The autocorrelation times in 
our simulations for iV = 20 and 40 are almost two times less than the ones for the same 
numbers of bosonic fields in Liischer's formulation measured in Ref . |Q , but our MC runs 
were rather short and that may tend to underestimating the algorithm's autocorrelation 
time. A more careful investigation of autocorrelation behaviour in Slavnov's algorithm 
requires further extensive tests on a more powerful computers. 

It is interesting that the systematic errors (9(6^) and 0{^^) from the one side and 
0{e~^^) from the other side for the measured quantities compensate each other. This 
fact can lead to the considerable improvement of results even for moderate values of 
as it was observed in the present research. 

The present performance of Slavnov's algorithm can be considerably improved by 
diminishing the systematic error from 0{h'^) to 0{h'^). This can be put into effect by 
adding to the action (|) the terms 

^S = Y. [«(6, AT, ii)x*{x)B\{x) + 5{h, N, fi)x*{x)x{x)] (26) 

X 

where the coefficients a{b, N, jj) and 5(6, A^, /i) can be computed theoretically (this work is 
under progress). Also even-odd preconditioning can be incorporated to reduce the number 
of points in auxiliary dimension A^. 
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